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Large biomolecular systems, whose function may involve thousands of atoms, cannot easily be 
addressed with parameter-free density functional theory (DFT) calculations. Until recently a central 
problem was that such systems possess an inherent sparseness, that is, they are formed from compo- 
nents that are mutually separated by low-electron-density regions where dispersive forces contribute 
significantly to the cohesion and behavior. The introduction of, for example, the van der Waals 
density functional (vdW-DF) method [PRL 92, 246401 (2004)] has addressed part of this sparse- 
matter system challenge. However, while a vdW-DF study is often as computationally efficient as a 
study performed in the generalized gradient approximation, the scope of large-sparse-matter DFT 
is still limited by computer time and memory. It is costly to self-consistently determine the electron 
wavefunctions and hence the kinetic-energy repulsion. In this paper we propose and evaluate an 
adaption of the Harris scheme [PRB 31, 1770 (1985)]. This is done to speed up non-selfconsistent 
vdW-DF studies of molecular-system interaction energies. Also, the Harris-type analysis establishes 
a formal link between dispersion-interaction effects on the effective potential for electron dynamics 
and the impact of including selfconsistency in vdW-DF calculations [PRB 76, 125112 (2007)]. 

PACS numbers: 31.15.E-,31.15.ae, 71.15.Nc 



I. INTRODUCTION 



Density Functional Theory (DFT) is considered one 
of the best and most reliable condensed-matter tools 
for non-empirical studies of molecular, surface, bulk and 
compound properties.^ Standard implementations, using 
the Local Density Approximation (LDA) or the General- 
ized Gradient Approximation (GGA) for the exchange- 
correlation energy, provide an accurate description of the 
binding in regions characterized by high electron density. 
An even more widespread DFT usage will follow from an 
ability to address soft- and sparse-matter systems, struc- 
tures that have internal voids or low-electron-density re- 
gions dominated by the van der Waals (vdW) forces, also 
called the London dispersion forcesP While neither LDA 
nor GGA capture the truly nonlocal correlation effects 
that underpin those forces,™^ the last decade has seen de- 
velopment of both vdW-extended DFT^HII] anc i f re g_ 
ular nonlocal exchange-correlation functionalsP^HIH xhe 
first class of methods are often atom centered and require 
use of a damping function or equivalent, while the second 
class of methods fits inside the regular DFT framework. 
Both types of sparse-matter DFT can describe, for ex- 
ample, vdW forces between molecules. 

The van der Waals density functional (vdW-DF) 
method 6 -* 7 * 13 * 14 * 19 !^ is a framework for approximat- 
ing the exchange-correlation energy E xc [n\. The 
method, summarized below, yi elds effi cient general- 
purpose sparse-matter functional d 12 * 13 * 19 ^ that are non- 
empirical. The method employs the Coulomb gauge 
(with Green function G = |r — r'| _1 ) and a scalar di- 
electric function e. In its most general formp^H the 
vdW-DF meth od is a reformulation of the adiabatic con- 
nection formuhPlU (ACF) and assumes that a plasmon- 



pole approximation^ 3 -* for e can satisfy 

J ^Tr{ln(VeVG)} = E xc [n] + E sel{ , (1) 

where u denotes the complex frequency and where E ae n 
is the internal Coulomb self-energy of each electron. This 
equation summarizes an, in principle, exact description 
of the (longitudinal) electrodynamics in the inhomoge- 
neous electron gas and reflects the use of a Dyson equa- 
tion for handling screening. Like the GGAs, th e vdW - 
DF method further uses physics-based constraints^^] t 
approximate the plasmon-pole response and thus defines 
the functional form of a nonlocal correlation term, E% 1 . 
In t he re cent explicit functional versions, called vdW- 
DFp2H2and vdW-DF2j±2l the nonlocal energy Ef[n] is 
expressed as a double integral over the density, weighted 
by a kernel. However, the plasmon basis still allows 
it to capture a collecti vity that reflects the broader 
density variation! 13 1 20 1 21 1 The vdW-DF method also in- 
volves picking a gradient-corrected exchange that re- 
flects prespecified cr iteria, e.g ., good all-round molecular- 
system p erform anc e 1 7 * 13 * 1 -^^* or improved bulk-system 
properties) 25 * 26 * The vdW-DF shares the plasmon-pole 
emphasis with its LDA and GGA relatives, and the func- 
tionals "vdW-DF#" have both seamless integration in 
the homogeneous limit and a build-in conservation of 
the exchange-correlation holeP^ The non-empirical de- 
sign suggests that the vdW-DFs can achieve a good trans- 
ferability across systems, length scales, charging states, 
and binding morphologies. 

The vdW-DF method has been and is being tested 
for many systems. It delivers a parameter- free atomic- 
scale characterization of the binding in complex sparse- 
matter systems. Selfconsistent (sc) vdW-DF calcula- 
tions can be used to calculate stress within periodic unit 
cells^and guide atomic optimization (relaxation). There 
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are performance tests for bulkj^HSSl layered^ 2 ' 34 1 35 1 
absorption pSMl molecule and atom adsorptionpS HUl su r _ 
face self-assemblyJ^2HHl an d molecular systems .H5156H1SI 
The development of effective algorithmsEHSSl has allowed 
sc vdW-DF calculations that have today little or no addi- 
tional costs over GGA calculations^ For very large sys- 
tems, typical of biomolecular-interaction problems, the 
set of real-space vdW-DF evaluation schemes^HsU effec- 
tively permits an order-TV scaling in the evaluation of 
[n] , as discussed for example in Ref . [3TJ 
Nevertheless, the feasibility of sparse-matter DFT cal- 
culations for large molecular matter is still limited by the 
huge number of atoms that are usually involved in inter- 
acting molecular systems. The problem is confounded by 
the fact that (especially in biochemically relevant prob- 
lems) we are often forced to address nonperiodic systems 
where size-conve rgenc e of the model unit-cell size be- 
comes important P^Hil Molecular t ransp ort studies,^ and 
many molecular-crystal problems] 31 * 67 ! further exemplify 
the general need to address vdW binding in big systems. 
The challenge for typical sparse-matter DFT is increas- 
ingly becoming that of computing the steric hindrance 
effects that are described by the DFT kinetic-energy re- 
pulsion. When we reach thousands of atoms, all first- 
principle and vdW-extended DFT calculations are sim- 
ply limited by memory requirements and computational 
costs of the wavefunction evaluation. 

In this paper we propose and test an adaption of the 
Harris scheme^^U t perform non-selfconsistent (nsc) 
vdW-DF calculations for sparse matter. The aim is to 
explore possibilities for reducing the computational cost 
of the wavefunction-evaluation bottleneck that could be 
impeding an even broader vdW-DF usage in, for exam- 
ple, biochemistry.^ The approach consists of using a 
superposition of the fragment (electron) densities, and 
is here termed sfd-vdW-DF. It can be seen as an al- 
ternative to the mo re cost ly nsc-vdW-DF or even sc- 
vdW-DF evaluations! 7 * 13 * 1 ^ It is just the regular vdW- 
DF Harris scheme if the fragment densities are based on 
vdW-DF calculations. However, the sfd framework also 
works without a vdW-DF implementation of the under- 
lying DFT code, and we reserve the term sfd-vdW-DF 
to cases where the fragment densities are obtained from 
GGA calculations. The method uses our code for real- 
space vdW-DF evaluation base d on the c harge densities 
from the underlying DFT code! 28 ' 31 ' 58 ' ^ 

We present formal analysis and a testing of the sfd- 
vdW-DF computational scheme for molecular systems. 
As an interesting aside, the analysis also identifies con- 
ditions for expecting significant vdW-induced changes 
in the bandstructure or electron dynamics for a given 
binding morphology. We separately test how well the 
sfd framework faithfully reproduces interaction effects 
that arise from the semi-local part of the vdW-DF func- 
tional. We test the performance of the sfd-vdW-DF 
scheme across the S22 benchmark setP^ and for systems 
with a varying degree of static polarizations. We ob- 
serve that while a vdW-DF Harris scheme (evaluated 



with vdW-DF fragment densities) is correct to second or- 
der in binding- induced density changes, Sn — n sc — n s fd, 
the proposed scheme is formally only correct to linear 
order in these density changes. However, we also find 
that the linear term Sn is weighted by the changes in the 
effective Kohn-Sham (KS) potential that result from the 
shift from the GGA to the vdW-DF functional. The sfd- 
vdW-DF scheme may thus be broadly applicable in the 
absence of large static dipoles. 

The Paper is organized as follows. In Section II we 
summarize the case for developing an accelerated (but 
approximate) vdW-DF description of biomolecular inter- 
action. In Section III we present some details of the LDA, 
GGA, and vdW-DF family of constraint-based density 
functionals to facilitate a formal analysis and our pro- 
posal for the sfd-vdW-DF scheme in Section IV. Section 
V provides a brief summary of computational (and im- 
plementation) details. Section VI presents the results of 
the performance testing that we have carried out for the 
proposed sfd-vdW-DF scheme, while Sec. VII contains a 
discussion. Finally, Section VIII contains summary and 
outlook. 



II. A CASE FOR FAST VDW-DF STUDIES OF 
BIOMOLECULAR INTERACTIONS 

In a related worlPSl we report a vdW-DF mapping 
of the vdW attraction in a DNA dimer (two peri- 
odic double- helix strands). That pilot study illustrates 
that general conclusions can be made for the important 
biomolecular-interaction problem from computing merely 
how E^ 1 depends on the interaction geometry. 

The nonempirical nature and the regular-density- 
functional basis (no external parameters) makes vdW- 
DF well suited to pursue investigations of biomolecular 
interactions. The full DNA interaction problem touches 
on two general challenges for refining a computational 
description of life processes. First, it illustrates the 
workings of molecular recognition (the matching by weak 
forces of the genes or just the packing of our genome in 
its environment). Second, it reminds us of the challenge 
of characterizing these effects in a solution that contains 
counter ions. Since the counter ions, per se, can be ex- 
pected to play a smaller role in the overall interfragment 
vdW attraction, it is natural to focus a sparse-matter 
DFT study on the charged biomolecular structures them- 
selves. However, for this approach to become meaningful, 
we must be able to also characterize the vdW bonding 
at various charging states. Being a regular (parameter 
free) nonlocal density functional, vdW-DF has an inher- 
ent advantage for computational studies of biomolecular 
interactions. 

The prospect of vdW-DF studies of biomolecular sys- 
tems is as promising in terms of computational cost (and, 
currently, as restrained) as it is for other sparse-matter 
DFT approaches. The evaluation of the nonlocal corre- 
lation, E~, causes no relevant slow down or bottlenecks 
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FIG. 1: Binding energy curve for the ammonia dimer, a 
member of the S22 benchmark set.^ The dark (light) solid 
curve indi cates the results for nonself consistent vdW-DFl^ 
(vdW-DFi^EH). The dashed dark (light) curve indicates the 
results of sfd-vdW-DFl (sfd-vdW-DF2). The cross identifies 
the binding energy and separation as obtained for the original, 
fully selfconsistent vdW-DF2 studyP 

thanks to the order- AT scaling in the real-space evalua- 
tion approach for large systems.^ In fact, our mapping 
of DNA attraction 7 ^ demonstrates that vdW-DF has, in 
practice, an efficiency similar to that of DFT-D for large 
biomolecular problems. There are for the DNA-dimer 
attraction problem^ no memory bottlenecks and essen- 
tially a linear scaling at least up to 1000+ cores for the 
E^ 1 evaluation. It is already today possible to bring that 
functional evaluation to about ten minutes wall time. 
While this evaluation of E™ y does cost more than adding 
the dispersion term in DFT-DpSl neither of these descrip- 
tions of the vdW attraction cause limitations: The com- 
puting expense is irrelevant when compared to the cost of 
converging the density in even one DNA (one periodically 
repeated coil of a DNA double helix) in DFT. 

With such a promise, it is frustrating that a full DFT 
study of the DNA dimer is today impossible without an 
allocation at a petafiop facility. There is a need to ac- 
celerate the DFT determination of large-system wave- 
functions and for sparse matter investigations in general. 
Here we explore an approach that focuses on accelerating 
the evaluation of the kinetic-energy repulsion, and which 
can, in principle, reduce the wavefunction-solution stage 
to just a single electroni c iterati on, as is possible with 
the regular Harris scheme. ^ 8 ' 70 ' 77 ' 

Figure [T] illustrates the feasibility of using the sfd-vdW- 
DF for accelerated vdW-DF studies of biomolecular in- 
teractions. It reports a comparison of binding in an am- 
monia dimer and shows that the sfd-vdW-DF2 descrip- 
tion is in excellent agreement with the binding predicted 
by sc-vdW-DF2 (indicated by a cross). 

Figure [2] summarizes our overall assessment, further 
detailed in Sec. VI. It testifies to a h igh degree of ro- 
bustness across the S22 benchmark setP^ and for other 
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FIG. 2: Organic-molecule assessment of the here-described 
(Harris-type) superposition of fragment densities (sfd) evalua- 
tion scheme, sfd-vdW-DF. Results for vdW-DFl (vdW-DF2) 
are compared in the top (bottom) panel. Each point repre- 
sents a set of binding energies for a pair of molecules in the 
S22 data set. The binding energies are for sfd-vdW-DF and 
sc-vdW-DF for vdW-DFl and vdW-DF2, and as calculated 
in quantum chemistry (QC) studies. The sc and QC results 
are listed in Ref. [19] QC results are originally from Ref. 1961 
All sfd and sc calculations are carried out at the sc-vdW-DFl 
or sc-vdW-DF2 binding distances.^ We note that the perfor- 
mance of sfd-vdW-DF is excellent for S22, both as compared 
to sc-vdW-DF and to QC calculations. 

molecular-interaction problems. The important observa- 
tion for developing a biomolecular computational strat- 
egy in DFT is that the sfd-vdW-DF scheme ca n be as re- 
liable as the often-used nsc-vdW-DF evaluation^- 3 ' 28 1 31 1 58 1 
even if it bypasses the need for a sc determination of the 
wavefunctions. The sfd framework could thus be one ap- 
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proach to speed up vdW-DF studies of large biomolecular 
interaction problems at a limited cost in accuracy. 

The (regular) Harris scheme works within a fixed func- 
tional choice and takes the following steps: First, self- 
consistent calculations of the densities 7ii=i,2,... for each 
individual building block second, construction of a 
density n s fd = J^. rii as a superposition of the build- 
ing blocks, and of the effective single-particle potential 
Kff.sfd(r) = Vrff[n s fd](r) that corresponds to n s fd; and 
third a nsc, that is, no-density update, calculation of 
the eigenvalues corresponding to the potential V^sfd (r) • 
These eigenvalues help provide an estimate of the kinetic 
energy term in the Harris scheme!^ The Harris scheme is 
traditionally pursued in a LDA or GGA framework start- 
ing, e.g., from sc GGA input densities. It is today often 
used when pursuing bandstructure calculations in DFT 
(and therefore provides very accurate wavefunctions for 
a given V e s, s {d(r)) ■ However, it is also possible to use it 
for its original purpose, namely for providing efficient but 
approximate DFT studies of interactions. 

Pursuing a regular vdW-DF Harris scheme represents 
one alternative for accelerating vdW-DF studies. We 
point out that the original study by Harris^ (working 
with LDA) shows that the scheme works reasonably well 
also for describing the formation of covalent bonds be- 
tween atoms in some molecules. It should be even better 
suited to describe the weaker vdW bonding. The over- 
all criteria for the applicability of a Harris-type scheme 
is that the fully sc density solution n sc should not dif- 
fer significantly from the input superposition density 
n. 8 f(i. It can work also when we consider systems, like 
biomolccules in solution, where the charging from the 
surrounding counter ions must be considered when pro- 
viding the input densities (and hence n s fd and V e s, s fd( T ))- 
A vdW-DF Harris scheme can be expected to work well 
for a study of supramolecular systems as long as the 
binding-induced charge relocations remain small. 

A regular Harris vdW-DF scheme is, however, not the 
focus here. This is because we do not yet have the abil- 
ity to both perform vdW-DF calculations and allow an 
externally-defined superposition of input density in the 
same code. We have chosen to build on the DacapcP^ 
code (for input densities and for the nsc kinetic-e nergy 
evaluation^) and on the vdW-DF postprocessing^^ 
that we have previously used extensively for vdW-DF 
studies. A benefit of introducing sfd-vdW-DF as a sup- 
plement to a regular vdW-DF Harris scheme is that we 
thus provide a computational framework that can take 
input densities from an arbitrary code (through an adap- 
tion of any o f the real-space strategies for evaluation 
vdW-DFEEUEI]) 



III. DENSITY FUNCTIONAL DESCRIPTIONS 
OF DENSE AND SPARSE MATERIALS 

We begin with a description of differences between the 
regular form of the nonlocal functional vdW-DF and the 



traditional local (LDA) and semilocal (GGA) DFT de- 
scriptions. This facilitates our subsequent discussion of 
the approximation that can allow an accelerated evalua- 
tion in (nsc) vdW-DF studies. 

A. Nonempirical approximations to the universal 
energy density functional, LDA, GGA, vdW-DF 

Both GGA and vdW-DF arc refinements of the ear- 
lier approach LDA that describes the universal exchange 
correlation functional 

*t° A M = E^ A [n]+E^ A [n] (2) 

*£,c A M = / d 3 rn(r)C(n(r)) (3) 

in terms of exchange and correlation energy densities 
e^ A (n(r)) that are just functions of the local density 
n(r). The behavior of the LDA was initially estab- 
lished by considering the self-energy shifts that result 
with a single-particle coupling to the collective plasmon 
excitations.^21221 The LDA description was later refined 
by considering Quantum Monte Carlo studies of the ho- 
mogeneous electron gas and its responsePsHSl 

The GGA adds a functional dependence on the local 
gradient through dimcnsionless measures of the gradient, 
s(r) = | Vn\/ (2kpn) and t(r) cx |Vn|/(nfc s ), where kp — 
(37r 2 n) 1//3 is the Fermi wave number, k s = (4kp /ttclo) 1 / 2 
denotes the Thomas-Fermi screening wave number, and 
ao = ti 2 /me 2 . The set of GGAs expresses the exchange- 
correlation energy 

Eg GA [n] ^|d 3 rn(r)[4(n(r), S (r))+ £ f(n(r),t(r))]; (4) 

we keep a subscript 'g' on the GGA energy densities to 
stress that one must pick a particular design choice, al- 
though history has pulled towards a few major choices.^ 
The inclusion of the dimensionless gradient permits 
a richer variation of functional forms. The set of 
constraint-based GGAs is among the most successful 
and these extend the plasm on picture of the LDA via a 
wavevector analysis! 23 * 82 * 83 ! while also emphasizing con- 
servation of the exchange-correlation hole. The develop- 
ment led to robust and very versatile GGA forms like the 
PBglH j n t ne constraint-based GGAs, the form of the 
exchange energy density 

4Wr), S (r)) = e^ A (n(r))F°( S (r)), (5) 

is given by an enhancement factor F^(s(r)) that has been 
chosen to satisfy a number of scaling laws, formal con- 
straints and guidelines JHsMH Formal analysis also guides 
the choice of the gradient-corrected correlation energy 
density eg (n(r), t(r)). It is important to note that the 
richer variation that Q supports is thus tempered by 
adherence to fundamental physics criteria and that the 
identification and use of these criteria has helped propel 
the GGA (and DFT) to a tremendous successP 
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The broad class of materials and systems that are char- 
acterized by sparseness does, however, require further re- 
finements beyond GGA. For example, organics, biomat- 
ter, and supramolecular systems are sparse in the sense 
that they have internal electron voids or other regions 
with a low electron distribution. Here the binding and 
function are dominated by the vdW forces that reflect an 
electrodynamic coupling and act across internal voids. 

The vdW-DF method goes beyond LDA and GGA 
by introducing a truly nonlocal correlation contribu- 
tion E^ 1 \n] that makes the electrodynamical coupling 



explicilfeli^i3lH v i a Eq. (|fj). The vdW-DF exchange- 
correlation energy is thus written 



E£*-™[n] = E£[ 



where 



E 



v,0t 



E. 



LDA r 



El 



(6) 



(7) 



denotes the semilocal part of the functional. W e use a 
superscript V to identify the vdW-DF version^- 12 * 13 * 19 ! 
and stress that these in general have different exchange 
components E%[n]. We note that the different vdW-DF 
versions also have different forms of the nonlocal correla- 
tion term but we have chosen not to make that explicit 
in our discussion. 

In the recent vdW-DF versions the nonlocal correla- 
tion term is expressed b y a se cond-order expansion in 
the plasmon-pole respons e 1 13 * 19 ! 



E 



ni r 



d 3 i 



d 3 r'n(r)0[ra](r,r')ra(r'). (8) 



This nonlocal correlation term still captures the broader 
density variation through a collectivity that the plasmon 
reflects and ([8| is designed so that the vdW-DF method 
avoids double counting with the terms captured in the 
local correlation. The vdW-DF# have a seamless inte- 
gration 



E:[n)=E^[n]+E?[n), 



(9) 



and thus bypass the need for using a damping function. 
As stressed in the introduction, the vdW-DF# are de- 
rived as an approximation to the ACF. The recent more 
explicit functionals also use the AC F to link the plasmon 
pole to an inner functiona l 13 * 14 * 1 9 * ^D that is also of the 
form (|7|). All functionals in the vdW-DF method build 
the nonlocal functional from the inside out, i.e., describe 
the electrodynamics of the dispersion forces by linking to 
response of the time-tested LDA/GGA plasmon descrip- 
tion. 

In this paper we use the vdW-DF method and work 
with both the vdW-DFl version (in which the revPBE 86 
GGA is used for £^ dW - DF [n]) a nd th e vdW-DF2 ver- 
sion (in which the re-fitted PW8tfM2U GGA is used for 
E x dw ~ BF [n)). We note that in additio n to the canonical 
Rutgers- Chalmers vdW-DF versions ^ 14 ' 19 ' 25 1 there are 
also variants! 88 * 89 ^ These variants fit the outer exchange 
functional £^ dw " DF [n] to a form that is fitted to specific 
data sets, for example the S22. 



B. The exchange-correlation potentials 

For a discussion of the nature of the Harris and 
nsc-vdW-DF schemes (below) it is important to also 
characterize differences in the corresponding exchange- 
correlation potentials 



/4cM(r) 



5n(r) ' 

^vdW-DF [ 

5n(r) 



(10) 
(11) 



We also use these potentials in a discussion of the error 
in the sfd-vdW-DF. Again we have used superscripts g 
and v to stress that for calculations we must pick specific 
versions of the GGA or of the vdW-DF. 

Ref . [14] provides a derivation of the effective-potential 
contribution A/x" 1 [n](r) that arises from taking func- 
tional derivatives of the nonlocal correlation term 
Ef. Fr om |6j) it follows that the vdW-DF exchange- 
correlation potential will also differ from a GGA 
exchange-correlation potential 



n) (r) 



A^[n](r) 



(12) 



by a semilocal potential term A//JJ> c °[n]. This semilocal 
potential change arises in part because vdW-DF sub- 
tracts off the gradient corrections to correlation. Also 
when discussing the difference from a given GGA 'g' 
of exchange-energy form E%. [n], the semilocal potential 
change must reflect the differences E v x [n] — E% [n] . 



C. Selfconsistent DFT calculations of the total 
energy 



The KS energy can be writterpS 
~1 



E KS [n]=T + J d 3 rn(r) 



<Wr) + K xt (r) 



+E xc [n]+E_ 
(13) 

Here K x t(r) is the external potential and </>n(r) the elec- 
trostatic potential at the given density n(r), 



3„/ n ( r ') 



(14) 



The first term Tq and the last term En of ( 13 ) express 



the kinetic energy of an effective single-particle wavefunc- 
tion problem and the internuclear repulsion term, respec- 
tively. 

Fully selfconsistent DFT calculations in the KS scheme 
proceed by finding the single-particle wavefunctional so- 
lution of an effective eigenvalue problem 



-^V 2 + V eS [n}{r) - e x \ ^> A (r) =0 ( I V) 
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defined by the density-dependent (and density-functional 
specific) effective potential 

KffN(r) = V ext (r) + <t> n {r) + ,u xc [n](r). (16) 

We use atomic units in all formal discussions of DFT cal- 
culations and the set of approximations. Selfconsistency 
is enforced by requiring that the resulting single-particle 
description of the electron density 



The semilocal functional component 



A 



(17) 



actually coincides with the density that specified the ef- 
fective single-particle potential (15). 



As an example of the total-internal energy DFT cal- 
culation that is thus made possible, we consider a two- 
fragment system with components separated by a dis- 
tance d. We here follow the presentation in Ref. [55] so 
as to simplify our subsequent discussion (Section IV). 
Summing up the set of occupied, single-particle energies 
e\, leads to an incorrect counting of the total electron- 
electron interaction energy. However, this complication 
is easily adjusted for, giving^ 

OCC „ s ^ 

E K s{d) = ^ eA ~y d3rn ( r ) j^scM +feW(r) 



+ E xc [n sc ] + E N (d). 



(18) 



All terms depend on the mutual fragment separation d 
(although we do not make that an explicit statement 
for all terms). We use £ ; KS W " DF ( d ) or E Ks( d ) to specify 
whether the sc DFT total energy result was pursued in a 
vdW-DF or a GGA choice, respectively. Below we focus 
the discussion on such fragment problems with mutual 
separation d. 



D. Non-selfconsistent vdW-DF calculations 

We present an overview of the standard nsc-vdW- 
DF evaluatiorPES w hich has a total-energ y va riation 
#nsc W-DF (d)- This energy variation is ofterPHi found 
to be a close approximation to the energy variation 
£vdW-DF( d ) = £vdW-DF( d ) found by fully sc v dW-DF 

calculations, Sec. III.C. 

The nsc-vdW-DF evaluations proceed for a given GGA 
'g' by first completing self-consistent GGA calculations of 
both the electron density variation n^ c (d) and total GGA 
internal energy E g (d). The GGA choice g is in practice 
often PBE and perhaps more seldom the revPBE version 
(that vdW-DFl uses for its exchange component but that 
is of no concern in this formal discussion) . We denote by 
E% c (d) the corresponding exchange and correlation en- 
ergy that are evaluated for nf c (d) . The nsc-vdW-DF cal- 
culations proceed by simply adjusting for the nonlocal 
correlation and for the differences in semi-local correla- 
tion terms 



E^ BA \ni 



+ El[nl c ]-E9 c {d) (20) 



not only extracts the gradient corrections to correla- 
tions but also implements a possible adjustment in the 
gradient-corrected exchange description (as necessary). 

The nsc-vdW-DF calculations were for some time the 
only manner for completing a vdW-DF study: It permit- 
ted us to include van der Waals interactions in a compu- 
tational efficient parameter-free single-density functional 
DFTpJHH] The approach can be motivated, in part by a 
Harris-type description (as substantiated further below) 
but the quantitative extent of the approximation could 
only be tested when efficient implemen tation s of the sc- 
vdW-DF method were made available! 14 * 61 ^ The subse- 
quent testing showed that nsc-vdW-DF often captures 
most of the binding of sc-vdW-DFP^ 



IV. HARRIS-TYPE EVALUATION SCHEMES 

This paper formally proposes a computational strategy 
that combines nsc-vdW-DF calculations (above) w ith a 
further approximation inspired by the Harris scheme^ 7 -^ 
and oth er ear lier suggestions of using frozen fragment 
densities! 71 * 72 ! The approximation can limit the computa- 
tional costs for molecular systems because it reduces the 
quality required for the input density in the nsc-vdW- 
DF evaluation. It comes with an accuracy cost, which as 
expected is found largest for systems with a static polar- 
ization, but it can provide a speed up. 



A. Variational nature of the Harris scheme 

The regular Harris scheme rests ultimately on the vari- 
ational character of the KS formulation of the total en- 
ergy for fully selfconsistent DFT evaluations. The KS 
energy form acquires a minimum at the correct ground- 
state density n sc , 

E KS [n sc + Sn}= E KS [n sc ] + C(5n) 2 ; C> 0. (21) 

The Harris scheme rewrites the KS formulation of the 
total energy so as to avoid the need for updating the 
electron density in an estimate of the interacting energy. 
As mentioned in the introduction, the Harris scheme does 
not provide nor does it work with the sc density n sc , but 
rather utilizes a superposition density 



«- s fd = > , n 



(22) 



E 



vdW-DF; 



(d) = ^(d)+^K c (d)]+A^ u K c (d)]. (19) 



defined from sc determinations of the electron densities 
n sc i for each of the fragments of the weakly interacting 
system. 
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In a GGA study, for example, we can formally express 
the Harris interaction estimate^ 



occ 



- | d 3 r < fd (r) |^„« fd (r) + MLK fd ](r) 
+ ^cK fd ]+^W- (23) 



Again, d is the distance between the fragments and 
is the electrostatic potential defined by n 9 {d . In the Harris 



sfd 



estimate (|23|), the eigenvalues e? are the single-particle 



energies calculated within the Harris "one-shot" (no den- 
sity update) GGA calculation for the frozen input su- 
perposition density n 9 {d using the effective (Harris-GGA) 
potential 

Cfd(r) = CKdKr) = y cxt (r)+^ fd (r)+^ c [< fd ](r). 

There does, of course, exist a corresponding expres- 
sion for a Harris approximation to vdW-DF calculations, 
in which case the input density would be n^ dVl : "' 



n vdW-DF 



sfd 



The central step in the Harris scheme is the assumption 
that the density change n sc — n s f d produces only a small 
change in the effective potential, 



AVeff(r) = 0„ SC O) - 0n sfd (r) 

+ A4>sc](r) - /4 c [n sfd ](r), 



(25) 



so that one can expand the difference in KS and Harris 
estimates for the single-particle energy sum 



occ 

E 



occ „ 

^e A , = / d 3 rn sc {r)AV cB {r) + 0(r 



■ n siA ) 2 . 

The linear term cancels out corresponding linear terms 
in the expansion of E^Ap\ and in the calculation of the 
electrostatic potential.^ 

The resulting single-shot (no density update) DFT es- 
timate is also variational 



larrisKfd] = ^Harris [™sc] + 0(n sc - 7J sfd ) , (27) 



However, unlike a regular vdW-DF Harris inter- 
built from the 
and not 



E^P F (d) h 



action estimate 

sc GGA result for the entire system 
the superposition of sc-vdW-DF fragment densities, 



-.vdW-DF _ 
'sfd 



-.vdW-DF 



The nsc-vdW-DF 



approximation ^ns^" 1 ^ can, however, still formally be 
seen as a further extension of the ideas that underpin 
the Harris estimate Eq. ( 23 ) . 



To establish a formal relation between sc and nsc vdW- 
DF calculations, we consider the differences in sc results 
that arise as we replace a GGA choice g with a vdW- 
DF choice for the exchange-correlation functional. We 
introduce 



6n si 



vdW-DF 



I 9 



vdW-DF - (j>ni c 

il sc sc 

i „vdW-DFr vdW-DFi 
"r H-xc l IL sc J 



(28) 



/4cK]> (29) 



to identify the changes resulting in the density and in the 
effectivepotential, respectively. As in the original Harris 
analysis,^ we can consider both 5n sc and AV SC small, 
and thus giving rise only to linear changes 



E< 



vdW-DF 



E4 



E ,■ 



vdW-DFi 



E xc [n a sc ] 



d 3 r< c dW - DF (r)AV/ sc (r)(30) 

d 3 rfe sc (rK™K c ](r). 

(31) 



Simply extending the analysis behind the Harris esti- 
mate therefore yields the formal relation 



E, 



vdW-DF 



(d) « ^ s (d) + ^ 1 K c ] + A^ c dW - DF '°Kc 



+ AE V ™- BF ^ + 0{8n sc f 



(32) 



= E. 



vdW-DF 



(d) + AE; d ™- DF ^ +0(Sn sc ) 2 

(33) 



with leading-order correction term 

A^vdW-DF^s = 



d 3 rfe sc (r){^ dW - DF [< dW - DF " 



(r)-/4 c [< c ](rX>4) 



but it is not, in general, an extremum! 68 ! 70 ! 91 -^ This fol- 
lows because there is no consistency between the Harris 
scheme input density n s fd, and the single-particle electron 
density fl7| ) that results with the Harris-scheme effective 
potential V r c ff[n sfd ](r). 



B. Nature of and error in non-selfconsistent 
vdW-DF calculations 

The non-selfconsistent vdW-DF total energy 
_E^ d c w " DF (d) is an approximation to the fully self- 
consistent vdW-DF result E^ W ~ DF (d) = E^~ DF (d). 



C. A sfd-vdW-DF scheme for accelerated 
calculations of molecular interactions 

We propose to pursue a Harris-type vdW-DF scheme 
that is based on the superposition n^ fd = nf + n 9 , + . . . 
of GGA fragment densities but which approximates the 
vdW-DF total energy by 



E. 



vdW-DF/'^ _ 
sfd 



(^ Harris W+^c 



sfdJ 



vdW-DF,0r 
xc L"sfdJ 

(35) 

Here, again, ^Hams(^) denotes the regular Harris esti- 
mate as described in the given GGA choice g. 



A formal analysis motivating the proposed vdW-DF 
approximation ( 35 1 is essentially already stated in Sec- 



tion III.B. We now consider slightly different density and 
effective-potential differences 



5n s fd 
AVrfa 



vdW-DF g 
n sid n s fd' 



,vdW-DF 
sfd 



sfd 



,,vdW-DFr vdW-DFi 
' Vxc [ n sfd J 



/4 C [ 



'sfdJ ' 



(36) 



(37) 



so that the g — > vdW-DF changes instead reflect the effects 
on the Harris-scheme single-particle cigcnenergics e^. 

The sfd-vdW-DF estimate can thus be expressed as an 
approximation to a regular vdW-DF Harris scheme 

= K f r DF (rf) + A< S ^ DF ^ 9 + 0{5n sc f 

where we now have a slightly different leading-order cor- 
rection term 



AE 



,vdW-DF<-g 
/x,sfd 



d 3 r5n sid (r) {/i 



vdW-DFr vdW-DF 



sl'd 



1«-Mf c [< fd ](#) 



V. COMPUTATIONAL DETAILS 

This paper compares the interaction energy curves 
obtained with nsc Harris-type calculations with those 
of DFT calculations for selected non-covalently bound 
molecular dimers. Both sc and nsc calculations with the 
PBE versiorPS of GGA are performed using the Dacapo 
software!^ This planewave DFT code was chosen because 
it is straightforward in Dacapo to set the electronic den- 
sities equal to the sum of molecular (frozen input) densi- 
ties n s fd = n,i + n 2 through an external manipulation in 
AS022I and thus to prepare the sfd calculations. 

The non-local correlation energy is evaluated in a post- 
processing procedure both for regular DFT and Harris- 
type vdW-DF calculations. For these calculations, we use 
an efficient in-house real-space code, further described in 
Ref. (3TJ A radius cutoff of 6 A is used for dense (full) 
sampling of the grid and a cutoff of 26 A is used for sparse 
(double-spaced) sampling of the grid. 

In the PBE calculations, relying on Vanderbilt ultra- 
soft pseudopotentials, we use plane-wave and density- 
sampling cutoffs of 500 eV. Thi s cutoff choice has been 
used in many similar calculation s] 30 * 31 ^ and gives a rela- 
tively dense sampling of the density grid used to evaluate 
the non-local correlation. As long as the reference calcu- 
lations have the same grid-sampling density, here secured 
by using the same size of the unit cell, the non-local cor- 
relation energy is typically converged to within about 
1 meV. 
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FIG. 3: A comparison between the regular Harris scheme and 
self-consistent (sc) Kohn-Sham calculations for four different 
molecular pair configurations identified in the upper panel, 
all systems investigated in the PBE version of GGA. These 
are (left to right, top to bottom) the HF dimer in parallel 
configuration, the HF-benzene pair with the H of HF point- 
ing towards the benzene center, the parallel benzene dimer, 
and the C60 dimer with the benzene rings facing each other. 
The middle panel shows the difference between the interac- 
tion energy in the sc and Harris PBE calculation. The lower 
panel shows (full curves) the interaction curves using the Har- 
ris functional and (dashed curves) the sc result. The abscissa 
label d denotes the separation between the closest atoms in 
separate molecules. The two curves involving the highly po- 
lar HF exhibit a non-negligible binding. The dotted vertical 
lines indicate the GGA-minimum. Among the investigated 
systems, the Harris estimate gives the largest overestimate 
(16%) for HF-benzene interaction curve. 



9 



VI. RESULTS: ASSESSING THE SFD-VDW-DF 
EVALUATION 

Four molecular pairs, depicted in the upper panel of 
Fig. [3| have been chosen for our comparison between the 
sfd scheme and regular DFT calculations. The first is 
a hydrogen fluoride (HF) dimer in parall el con figuration. 
This configuration is not the optimal onefSES but is here 
chosen as a representative for systems with large dipolc- 
dipole interactions. The second is a molecular configura- 
tion where the hydrogen of HF points towards the center 
of a benzene molecule. Thus, one molecule has zero and 
the other a large dipole moment in vacuum. The third 
system is a benzene dimer in parallel sandwich configu- 
ration. The binding in this system is dominated by vdW 
(also called London dispersion) forces. The interaction in 
this system is representative of dilute sparse matter sys- 
tem, like a gas. The fourth system, a dimer of C60 with 
hexagonal rings facing each other, is also one where the 
binding is dominated by vdW forces. But because of the 
large size of C60, this attraction is much stronger than 
for the benzene dimer. This system is therefore more rep- 
resentative of compact molecular complexes that arise in 
bulk sparse matter. 



A. Regular Harris scheme for GGA-PBE 
calculations 

We first describe and illustrate the Harris scheme as 
it is used for GGA calculations. Obviously, the use of a 
GGA will not generally succeed in reproducing structural 
properties of typical sparse, weakly interacting molecular 
systems. Nevertheless, it is instructive to illustrate that 
the GGA Harris scheme is still generally able of faith- 
fully reproducing the sc GGA calculations, including the 
sparse-matter GGA limitations. 

Figure [3] compares the interaction curves for the four 
different molecular pairs as obtained with DFT and 
the Harris scheme using the PBE version of the GGA 
exchange-correlation functional. Only the HF dimer and 
the HF-benzene pair show an appreciable binding of re- 
spectively 173 and 180 meV using the Harris scheme for 
PBE and 158 and 156 meV using sc PBE calculations. 
For the parallel HF dimer system, which is dominated by 
dipole-dipole interactions, the Harris calculation overes- 
timates the binding energy by 9% compared to regular 
GGA DFT calculations. For the HF-benzene system, 
where one of the molecules is highly polar and the other 
is not, the scheme overestimates the binding energy by 
17%. 

The discrepancies between the two methods can be 
understood from the significant dipole moment induced 
by the binding. At optimal separation (in the selected 
configurations), a dipole of 0.12 eA is induced for the HF 
dimer, while one of 0.15 eA is induced for the HF-benzene 
pair. These induced dipole moments are comparable to 
the dipole moment of the HF molecule itself (0.39 eA). 




FIG. 4: Comparison of the sfd-vdW-DF scheme for a benzene 
dimer and results obtained in nsc-vdW-DF evaluation. The 
top panel compares the differences between these two proce- 
dures as they arise for the total energy of vdW-DFl (black 
curve with star dots) and vdW-DF2 (cyan/light curves with 
star dots); and the difference in the nonlocal correlation (two 
upper curves). The bottom panel compares the sfd-vdW-DFl 
(black star-dotted curves) and sfd-vdW-DF2 (cyan/light star- 
dotted curves) against nsc-vdW-DF results (dashed curves). 
The somewhat larger shift (for most separations) in the non- 
local correlation energy when using the sfd scheme compared 
to the total vdW-DF energy, in particular for vdW-DF2, in- 
dicates a partial cancellation of the density sensitivity of the 
nonlocal and semilocal components. 



It is clear that molecular pairs involving one or more HF 
molecule(s) serve as tough tests for the feasibility of the 
Harris functional scheme. 

For systems dominated by the vdW forces, the discrep- 
ancies are difficult to assess without including the effect 
of non-local correlation. We will therefore make this as- 
sessment in the next subsection. 
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B. Systems dominated by vdW attraction 

The benzene dimer is a typical organic system bound 
by vdW forces. Since this system is weakly bonded, we 
can expect charge transfer to be small and thus the vdW- 
DF Harris scheme, and more generally the sfd framework, 
to be well suited to describe the system. 

Figure [4] shows the comparison between the sfd and nsc 
vdW-DF calculations for the benzene dimer. The differ- 
ence between the dashed and the full curves in the lower 
panel is barely distinguishable. At binding separation 
the sfd result is 2% below the nsc result for vdW-DFl 
and merely 0.4% for vdW-DF2. 

The upper panel of Fig. [4] shows that the non-local 
energy is somewhat affected by using the frozen density 
n s fd = m + 7i2 in place of the one determined with a full 
GGA calculation, n nsc . It also reveals that there is some 
error cancellation between these shifts and the combined 
shifts in the other terms: the shifts obtained with the 
sfd scheme overestimate the non-local interaction energy, 
while the magnitude of the binding energy is underesti- 
mated. This trend is opposite to that exhibited for all 
four systems in the proper GGA Harris scheme, shown 
in Fig. [3] For vdW-DF2 this error cancellation is close 
to exact in a fairly wide region around the binding sep- 
aration. For shorter separation between the molecules, 
corresponding to a larger density overlap, the absolute 
difference between schemes increases, as does the magni- 
tude of the repulsive wall between molecules. 

Figure [5] compares the two methods for a C60 dimer in 
the same fashion as for the benzene dimer. In this case 
the sfd-vdW-DFl underestimates the binding energy by 
as little as 0.2%, while sfd-vdW-DF2 is spot on (within 
about 0.05 meV). This striking coincidence (arising from 
error cancellation) is likely fortuitous since the results 
are similar, but not this similar, in other regions of the 
interaction curve. For the C60 interaction curve, the sfd- 
vdW-DF calculations in some regions overestimate and 
in other regions underestimate the interaction energy. 

The benzene and C60 dimer calculations indicate that 
the sfd scheme is an appropriate method to accelerate the 
evaluation of interaction energies in systems dominated 
by vdW interactions. 



C. System with large induced charge: HF 
interacting with benzene 

In the HF-benzene system the vdW forces contribute 
to the binding alongside electrostatic effects. The nsc- 
vdW-DF2 predicts a binding energy of 174 meV com- 
pared to that of 155 meV with the sc PBE calculations 
(in Fig. [3k. 

Figure 6] shows the results of the interaction curves 
obtained with sfd and nsc vdW-DF calculations. For 
this system, we also find that the vdW-DF2 calculation 
produces a larger binding energy than the vdW-DFl, 
which is opposite to the case for the benzene and for 
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nl.vdW-DFl 



nl.vdW-DFl 



m>dW-DF2 n,vdW-DF2 



nl.vdW-DF2 
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FIG. 5: The sfd-vdW-DF description of binding in the C60 
dimer, legends and details as described in Fig. [4] 



the C60 dimer. This switching of order is related to 
the fact that vdW-DF2 in general has a less repulsive 
exchange accounPS and a less attractive non-local cor- 
relation account Since the smaller size of this system 
decreases the magnitude of the non-local correlation, it 
shifts the balance between the repulsive and attractive 
terms. 

The difference between the sfd and nsc vdW-DF calcu- 
lations increases to as much as 30% for vdW-DF2, com- 
pared to 17% for the sfd and sc PBE calculations. The 
discrepancy is somewhat smaller for vdW-DFl. Note 
that this discrepancy arises mostly from the shift in the 
AK^? W ~ DF '° term and not from the non-local correlation, 
which contributes with 4 meV in the opposite direction 
of the total shift of —60 meV. 

Our results indicate that the increased inaccuracy of 
the sfd scheme for polar systems arises primarily from 
short-ranged effects. Thus, the inaccuracy may be re- 
duced by starting from densities generated with revPBE 
exchange for sfd-vdW-DFl and in the same vein PW86r 
for sfd-vdW-DF2. 
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FIG. 6: The sfd-vdW-DF description of binding between a 
HF molecule and benzene. Legends and details as in the 
lower panel of Fig. [4] The sizable discrepancy between the sfd 
and sc calculation in the PBE calculations, shown in Fig. [3] 
carries through to the case of (sfd and nsc) vdW-DF calcula- 
tions. The non-local correlation has only a tiny effect on the 
discrepancy. 



D. The nsc-vdW-DF approximation 



We note that the leading-order difference (34) be- 



tween sc- and nsc-vdW-DF total-energy results is nomi- 
nally linear in the density change, Sn sc — n^ w ~ BF (r) — 
7i| c (r). The nsc-vdW-DF calculations — while often very 
successful — need therefore not always be as robust as a 
regular Harris vdW-DF scheme would be. 

On the other hand, the regular nsc-vdW-DF approaclP 
does have a mechanism for including some of the electron 
density rearrangement that arises from Pauli exclusion or 
from the formation of more traditional types of bonding 
(those that a GGA does capture). 

Fig. [6] shows that keeping such charge adjustments can 
be important for systems where at least one fragment has 
a large static polarization. Generally, we expect the nsc- 
vdW-DF approach to be more accurate except in cases 
with very weak intermolecular interactions. Of course, 
the only way to resolve the diff erence w ould be to perform 
a fully sc vdW-DF calculationPSUHl The sc-vdW-DF is 
now becoming standard procedure for medium to large 
systems (system sizes approaching a thousand atoms). 
However, a performance testing comparing £^ C W-DF and 



T^vdW-DF 
^sfd 



against E s 



vdW-DF 



= E 



vdW-DF 
KS 



is for technical 



reasons beyond the scope of this paper (as discussed in 
Section II). 



E. An organic-molecular testing of sfd-vdW-DF 
performance 

Figs. [T] and [2] and Table [I] presented a summary of the 
further assessment we have performed of the accuracy of 



the sfd-vdW-DF scheme for the S22 benchmark suited 
Here we provide some additional details. 

Table H] presents the calculated numbers of our compar- 
ison of sc-, nsc-, and sfd-vdW-DFl and -vdW-DF2 results 
for interaction energies in the S22 set of molecular dimers. 
These interaction energies are all evaluated at the bind- 
ing distance (identified in Ref. 19) that minimizes re- 
spectively the sc-vdW-DFl and sc-vdW-DF2 interaction- 
energy variation. The quantum chemistry computations 
are from Ref. |96l For an actual S22-benchmarking of 
various sparse-matter DFT methods one should provide 
a full binding energy curve for each of the computational 
approaches. Here, our purpose is merely to complement 
our analysis based on binding curves of illustrative spe- 
cial cases with statistics for the S22 set of dimers that are 
seen as typical of organic- molecular interaction problems. 

We note that system 1-7 can be labeled hydrogen- 
bonding dominated, while 8-15 can be labeled dispersion 
dominated, and the remainder mixed. By studying the 
table, it becomes clear that sfd tends to compare well 
with nsc results for dispersion-dominated systems, while 
the biggest discrepancies arise among the systems dom- 
inated by hydrogen bonds. This observation agrees well 
with our analysis based on Figs. |3j|6| 

Figure [2] conveys an overview and feeling for the qual- 
ity of the sfd calculations compared to nsc and sc cal- 
culations. Together, the figure and table show that nsc 
and sc calculations are very similar. This is reassuring 
considering the fact that they are also based on different 
codes. As earlier discussed the sfd results compare well 
to the nsc results. Further, the inaccuracy introduced is 
overall smaller than the difference between sc-vdW-DFl 
and QC results, while the inaccuracy is about equal to 
the difference between QC and vdW-DF2 results. vdW- 
DF2 has better performance for the S22 data set than 
vdW-DFl. The inaccuracy introduced by using sfd does 
not necessarily make results compare worse to QC. 



VII. DISCUSSION 

A. On vdW-bonding effects on electron dynamics 

We begin with an interesting aside, noting the impli- 
cations of our formal analysis on the expected error in 
nsc-vdW-DF and in sfd-vdW-DF, Eqs. pil) and p9|. 



It is known that the inclusion of nonlocal correlation 
and v dW fo rces often gives rise to indirect bandstructure 
effect because the vdW binding changes the mor- 
phology and hence the local environment for the elec- 
tron dynamics. However, it is also interesting to identify 
conditions where one can also expect direct vdW band- 
structure effects, that is, electron-dynamics changes that 
arise when — for given structure — the effective potential 
is changed from a GGA to the vdW-DF functional form. 

The error estimate ( 34 ) allows us to identify conditions 



for expecting the vdW-bonding to affect the bandstruc- 
ture and, more generally, the electron dynamics! 33 * 97 * A 
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TABLE I: Interaction energies for pairs of small molecules from the S22 dataset. Quantum chemistry (QC) results from Ref. 
1961 the selfconsistent (sc) vdW-DFl and vdW-DF2 results are from Ref. 1191 All energies in meV/dimer. 



II 

# 


Duplex 


pvdW-DFl 


pvdW-DFl 


^vdW-DFl 


^sfd 


p vdW-DF2 

-^nsc 


c,vdW-DF2 


QC 


1 


Ammonia dimer 


118 


111 


115 


136 


126 


134 


137 


2 


Water dimer 


186 


171 


185 


220 


201 


218 


218 


3 


Formic acid dimer 


736 


680 


690 


817 


745 


766 


815 


4 


Formamide dimer 


610 


560 


587 


684 


619 


655 


699 


5 


Uracil dimer 


807 


749 


767 


876 


805 


832 


897 


6 


2-pyridoxine - 2-aminopyridine 


671 


623 


639 


728 


661 


687 


737 


7 


Adenine - thymine 


652 


619 


609 


716 


661 


660 


726 


8 


Methane dimer 


38 


39 


36 


29 


30 


30 


23 


9 


Ethene dimer 


70 


67 


64 


65 


61 


65 


65 


10 


Benzene - methane 


72 


70 


68 


64 


62 


63 


63 


11 


Benzene dimer (slip-parallel) 


136 


141 


136 


120 


124 


123 


114 


12 


Pyrazine dimer 


189 


188 


185 


178 


177 


177 


182 


13 


Uracil dimer (stacked) 


414 


396 


403 


412 


391 


402 


422 


14 


Indole - benzene (stacked) 


231 


232 


206 


234 


199 


197 


199 


15 


Adenine - thymine (stacked) 


466 


456 


461 


467 


457 


466 


506 


16 


Ethene - ethine 


74 


68 


69 


74 


67 


70 


65 


17 


Benzene - water 


142 


125 


124 


148 


126 


129 


143 


18 


Benzene - ammonia 


104 


97 


94 


101 


91 


92 


101 


19 


Benzene - HCN 


194 


162 


166 


198 


159 


170 


197 


20 


Benzene dimer (T-shape) 


120 


115 


113 


108 


101 


105 


118 


21 


Indole - benzene (T-shape) 


240 


210 


214 


230 


192 


206 


243 


22 


Phenol dimer 


267 


252 


254 


290 


270 


279 


307 



clear difference in sc and nsc vdW-DF total energies is 
required in order for such direct vdW-DF bandstructure 
effects to emerge. 

Meanwhile, there now exists a significant experience 
with using vdW-DF for sparse matter systems and the 
calculations have shown that there is, in practice, often 
only limited differences in nsc-vdW-DF and sc-vdW-DF 
calculations ISIMI follows that one must, in general, ex- 
pect that direct bandstructure effects typically are small. 



B. Overall assessment of the sfd-vdW-DF 
calculation scheme 

The speed up gained when using the sfd calculations 
with Dacapo are substantial yet somewhat modest. 
Computational costs are reduced by 40% and 55% for 
the benzene and C60 dimer respectively, when using stan- 
dard cutoffs with a minimal number of bands. Consid- 
ering that this software usually requires about 20 elec- 
tronic iterations to converge at these system sizes (but 
more for large systems), this gain is less than one might 
anticipate.^ 

However, we should consider that standard software 
like Dacapo (that we here use) has been subjected to 
intense efforts to optimize its ability to simultaneously 
solve the problem of charge relaxation and determination 



of the KS eigenvalues. What formally constitutes Harris 
calculations in that code are today primarily used to ob- 
tain accurate values for the KS eigenvalues.^ We do not 
desire such enhanced accuracy for an actual sfd-vdW-DF 
study. 

The fact that the here-proposed sfd-vdW-DF is still 
faster than nsc-vdW-DF (Sec. VI) is therefore promis- 
ing. Furthermore, since the performance, documented 
here, is excellent for many molecular systems, there is 
room for more compromise on accuracy. We believe that 
the present results motivate the approach to be further 
evaluated in forthcoming studies.^ 



C. Towards as fast biomolecular mapping of vdW 
interactions in biomolecular systems 

We are ultimately interested in vdW-DF computa- 
tional studies for large-scale interaction problems where 
there are relevant speed ups to be gained by not seek- 
ing the fully sc density (as described cither in a full 
GGA or in a full vdW-DF study). The acceleration must 
come from minimizing the cost of DFT calculations of 
the steric-hindrance or kinetic-energy repulsion effects 
(as this is the large-system bottleneck). 

We note that a parallel study, Ref. [751 pursues a 
closely related computational strategy and begins a first- 
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principle DFT mapping of the morphology variation on 
the vdW attraction in a DNA dimer by first focusing on 
an efficient evaluation of [«f fd ] ■ The DNA dimer sys- 
tem is there taken as an example of a typical large-scale 
biomolecular interaction problem. 

The vdW-DF exploration that is proposed in Ref. [76] 
bypasses the need for performing the expensive compu- 
tations of the kinetic-energy repulsion term (which, in 
any case, is not relevant outside the binding regime that 
extends to a nearest-atoms separation of about 4.0 A, 
Ref.HU). The first-principle vdW-DF survey of DNA at- 
traction thus achieves a dramatic speed up but can still 
(for relevant large-molecular system interaction geome- 
tries) be supplemented by adding the other components 
of the Harris-type sfd-vdW-DF scheme ( 35 ) . 



VIII. SUMMARY AND OUTLOOK 

To accelerate large-scale vdW-DF characterization of 
biomolecular systems it see ms useful to adapt the ideas 
of the Harris schem e 68 ! 70 " 72 ! as is indicated and ex- 
plored here. A Harris- type approximation which works 
reasonably well for describing the kinetic-energy effects 
of forming covalent bonds of atoms in some molecules 
should have a good chance of describing the simpler 
kinetic-energy repulsion (steric hindrance) of molecules 
in supramolecular systems. 

Here, we have put this expectation to the test. Our re- 
sults indicate that this scheme is promising for describing 
supramolecular systems bonded primarily by vdW forces. 
However, if one or more fragments are highly polar, this 
comes at the cost of accuracy. 

This paper is also supplemented by a related 
publication 76 which presents a vdW-DF study that maps 



out the nonlocal correlation of large biopolymers within 
the presented sfd-vdW-DF scheme. 

The pair of papers suggest a possible computational 
strategy for the study of binding in large supramolec- 
ular systems. The suggestion is to begin the structure 
and interaction-morphology search by essentially cost- 
free evaluations of the E^ 1 variation. That E~ step 
is available simply from relatively cheap calculations of 
fragment electron densities. One can in turn search for 
relevant binding motifs, given the linking to the here- 
proposed and tested sfd-vdW-DF scheme. This strategy 
eventually leads to a complete vdW-DF estimate of the 
variation in total interaction energy. The strategy exists 
as an alternative to implementing a real-space vdW-DF 
version in a code that realizes genuine ov&er-N scaling 
for large systems. 

Having established the promise of the sfd scheme for 
systems bound by vdW forces, the next step would be to 
investigate if the scheme can be further accelerated, in 
particular for large supramolecular systems. To this end, 
implementations for a full sfd-vdW-DF scheme for other 
DFT codes (beyond Dacapo) are being tested. 
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